1 Overview

This document contains all code necessary to replicate the analyses for the above article. Any questions about this project can be directed to Kaysee Arrowsmith at .

2 Setup

2.1 Load Packages

library(knitr)
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
library(tidyverse)
theme_set(theme_classic())
library(vegan)
library(ecodist)
library(mclogit)
library(glmmTMB)
library(DHARMa)
library(RColorBrewer)
library(geosphere)
library(Matrix)
library(lme4)
library(lmerTest)
library(ggeffects)

2.2 Matrix Regression Code

MRM.base <- function (formula = formula(data), data, nperm = 1000, method = "linear", mrank = FALSE, permute.control = "permute.control") {
  m <- match.call(expand.dots = FALSE)
  m2 <- match(c("formula", "data"), names(m), nomatch = 0)
  m <- m[c(1, m2)]
  m[[1]] <- as.name("model.frame")
  m <- eval(m, parent.frame())
  m <- as.matrix(m)
  n <- (1 + sqrt(1 + 8 * nrow(m)))/2
  if (abs(n - round(n)) > 1e-07) 
    stop("Matrix not square.\n")
  n <- round(n)
  if (ncol(m) < 2) 
    stop("Not enough data. \n")
  if (method == "linear") {
    X <- m[, 2:ncol(m), drop = FALSE]
      Y <- m[, 1, drop = FALSE]
      colnames(Y) <- "Y"
      newdata <- data.frame(Y = Y, X)
      fit1 <- lm(Y ~ ., data = newdata)
      b <- coefficients(fit1)
      sig <- summary(fit1)$sigma
      sig.df <- summary(fit1)$df[2]
      b.pval <- NA
      dev.pval <- NA
      if (nperm > 0) {
        b.all <- matrix(NA, nrow = nperm, ncol = length(b))
        b.all[1, ] <- b
        sig.all <- rep(NA, nperm)
        sig.all[1] <- sig
        for (i in 2:nperm) {
          newSample <- shuffle(n, control = permute.control)
          newY <- full(Y)
          newY <- newY[newSample, newSample]
          newY <- lower(newY)
          newdata <- data.frame(Y = newY, X = X)
          newfit <- lm(Y ~ ., data = newdata)
          b.all[i, ] <- coefficients(newfit)
          sig.all[i] <- summary(newfit)$sigma
        }
        b.pval <- apply(b.all, 2, function(x) length(x[abs(x) >= 
          abs(x[1])])/nperm)
        R2.pval <- length(sig.all[sig.all >= sig.all[1]])/nperm
      }
      results <- list(coef = cbind(b, pval = b.pval), 
        R2 = c(R2 = sig, R2.df = sig.df, 
          R2.pval = R2.pval))
  }
  else {
    if (method == "logistic") {
      X <- m[, 2:ncol(m), drop = FALSE]
      Y <- m[, 1, drop = FALSE]
      colnames(Y) <- "Y"
      newdata <- data.frame(Y = Y, X)
      fit1 <- glm(Y ~ ., data = newdata, family = binomial(link = "logit"))
      b <- coefficients(fit1)
      dev <- summary(fit1)$deviance
      dev.df <- summary(fit1)$df.residual
      b.pval <- NA
      dev.pval <- NA
      if (nperm > 0) {
        b.all <- matrix(NA, nrow = nperm, ncol = length(b))
        b.all[1, ] <- b
        dev.all <- rep(NA, nperm)
        dev.all[1] <- dev
        for (i in 2:nperm) {
          newSample <- shuffle(n, control = permute.control)
          newY <- full(Y)
          newY <- newY[newSample, newSample]
          newY <- lower(newY)
          newdata <- data.frame(Y = newY, X = X)
          newfit <- glm(Y ~ ., data = newdata, family = binomial(link = "logit"))
          b.all[i, ] <- coefficients(newfit)
          dev.all[i] <- summary(newfit)$deviance
        }
        b.pval <- apply(b.all, 2, function(x) length(x[abs(x) >= 
          abs(x[1])])/nperm)
        dev.pval <- length(dev.all[dev.all >= dev.all[1]])/nperm
      }
      results <- list(coef = cbind(b, pval = b.pval), 
        dev = c(resid.dev = dev, resid.df = dev.df, 
          dev.pval = dev.pval))
    }
    else {
      stop("method must be 'linear' or 'logistic'\n")
    }
    }
  results
}

2.3 Load Data

dat <- read.csv("data/tempint-networks.csv", stringsAsFactors = F) %>%
  mutate(date = as.Date(date, "%F"),
         year = format(date, "%Y"))
plants <- read.csv("data/tempint-plants.csv", stringsAsFactors = F) %>%
  mutate(date = as.Date(date, "%F"),
         year = format(date, "%Y"))
samptemps <- read.csv("data/tempint-samptemps.csv", stringsAsFactors = F) %>%
  mutate(date = as.Date(date, "%F"),
         year = format(date, "%Y")) %>%
  filter(!is.na(round),
         site %in% dat$site)
siteinfo <- read.csv("data/tempint-siteinfo.csv", stringsAsFactors = F)

3 Data Prep

3.1 Possible Interactions

Calculate every possible interaction that could have occurred based on the plants and pollinators observed at that date/site/round.

# Family
poss.family <- dat %>%
  filter(insect_family != "") %>%
  group_by(date, site, round, insect_family) %>%
  tally() %>%
  full_join(plants, by = c("date", "site"), multiple = "all") %>%
  dplyr::select(date, site, round, insect_family, plant) %>%
  mutate(concat = as.factor(paste(date, site, round, sep = ".")),
         date.site = as.factor(paste(date, site, sep = ".")))
  
# Genus
poss.genus <- dat %>%
  filter(insect_genus != "") %>%
  group_by(date, site, round, insect_genus) %>%
  tally() %>%
  full_join(plants, by = c("date", "site"), multiple = "all") %>%
  dplyr::select(date, site, round, insect_genus, plant) %>%
    mutate(concat = as.factor(paste(date, site, round, sep = ".")),
         date.site = as.factor(paste(date, site, sep = ".")))

# Functional Group
poss.func <- dat %>%
  filter(gross_ID != "Unknown") %>%
  mutate_if(grepl("Bombus*", .), ~replace(., grepl("Bombus*", .), "Bombus")) %>%
  group_by(date, site, round, gross_ID) %>%
  tally() %>%
  full_join(plants, by = c("date", "site"), multiple = "all") %>%
  dplyr::select(date, site, round, gross_ID, plant) %>%
    mutate(concat = as.factor(paste(date, site, round, sep = ".")),
         date.site = as.factor(paste(date, site, sep = ".")))

3.2 Realized Interactions

Which of the possible interactions (above) actually occurred.

# Family
occ.family <- dat %>%
  filter(insect_family != "") %>%
  group_by(date, site, round, insect_family, plant) %>%
  tally(name = "int") %>%
  right_join(poss.family, by = c("date", "site", "round", "insect_family", "plant")) %>%
  mutate(int = ifelse(is.na(int), 0, int))

# Genus
occ.genus <- dat %>%
  filter(insect_genus != "") %>%
  group_by(date, site, round, insect_genus, plant) %>%
  tally(name = "int") %>%
  right_join(poss.genus, by = c("date", "site", "round", "insect_genus", "plant")) %>%
  mutate(int = ifelse(is.na(int), 0, int))

# Functional Group
occ.func <- dat %>%
  filter(gross_ID != "Unknown") %>%
  mutate_if(grepl("Bombus*", .), ~replace(., grepl("Bombus*", .), "Bombus")) %>%
  group_by(date, site, round, gross_ID, plant) %>%
  tally(name = "int") %>%
  right_join(poss.func, by = c("date", "site", "round", "gross_ID", "plant")) %>%
  mutate(int = ifelse(is.na(int), 0, int))

3.3 Join Temperature

Attach the temperature for each sampling round. We also calculate some relative abundances here that will be useful for plotting our results later on.

# Family
final.family <- occ.family %>%
  left_join(samptemps, by = c("date", "site", "round")) %>%
  mutate(int.ID = paste(insect_family, plant, sep = ".")) %>%
  dplyr::select(date, site, round, date.site, concat, insect_family, plant, int.ID, int, avg.temp) #, latitude, longitude)

# Genus
final.genus <- occ.genus %>%
  left_join(samptemps, by = c("date", "site", "round")) %>%
  mutate(int.ID = paste(insect_genus, plant, sep = "."),
         year = format(date, "%Y")) %>%
         # midtime = as.numeric(strptime(starttime, format = "%H:%M")) + as.numeric(strptime(endtime, format = "%H:%M"))/2,
         # doy = format(date, "%j")) %>%
  dplyr::select(date, year, site, round, date.site, concat, insect_genus, plant, int.ID, int, avg.temp)#, latitude, longitude, midtime, doy)

final.func <- occ.func %>%
  left_join(samptemps, by = c("date", "site", "round")) %>%
  mutate(int.ID = paste(gross_ID, plant, sep = "."),
         year = format(date, "%Y")) %>%
         # midtime = as.numeric(strptime(starttime, format = "%H:%M")) + as.numeric(strptime(endtime, format = "%H:%M"))/2,
         # doy = format(date, "%j")) %>%
  dplyr::select(date, year, site, round, date.site, concat, gross_ID, plant, int.ID, int, avg.temp)#, latitude, longitude, midtime, doy)

4 Analysis

4.1 Relationship between Temperature and Round

This checks whether temperature varied with round. A priori, we would assume that temperature varies with date, but that should not matter in our analysis since our choice sets are always within a single date. However, we include date here too in order to control for that effect.

round.mod <- lmer(avg.temp ~ date + as.factor(round) + (1|site), data = samptemps)
summary(round.mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg.temp ~ date + as.factor(round) + (1 | site)
##    Data: samptemps
## 
## REML criterion at convergence: 1156.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.66524 -0.58990 -0.04593  0.65944  2.33962 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept)  2.275   1.508   
##  Residual             16.298   4.037   
## Number of obs: 203, groups:  site, 6
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       101.689518  31.697145 197.176783   3.208  0.00156 ** 
## date               -0.004049   0.001663 197.414224  -2.436  0.01576 *  
## as.factor(round)2   2.597828   0.680227 194.334156   3.819  0.00018 ***
## as.factor(round)3   3.253508   0.704645 194.432966   4.617 7.06e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) date   as.()2
## date        -1.000              
## as.fctr(r)2 -0.015  0.005       
## as.fctr(r)3 -0.022  0.012  0.487
anova(round.mod)
## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## date              96.68  96.684     1 197.41  5.9321   0.01576 *  
## as.factor(round) 400.27 200.134     2 194.37 12.2795 9.508e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
samptemps.BC21 <- samptemps %>% filter(site == "Brush Creek Hill", year == "2021")
test.bc21 <- aov(avg.temp ~ as.factor(round), data = samptemps.BC21)
summary(test.bc21)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2   6.32   3.158    0.18  0.837
## Residuals        14 245.65  17.547
samptemps.BC22 <- samptemps %>% filter(site == "Brush Creek Hill", year == "2022")
test.bc22 <- aov(avg.temp ~ as.factor(round), data = samptemps.BC22)
summary(test.bc22)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2   61.9   30.93   1.184  0.333
## Residuals        15  391.7   26.11
samptemps.KP21 <- samptemps %>% filter(site == "Kebler Pass", year == "2021")
test.kp21 <- aov(avg.temp ~ as.factor(round), data = samptemps.KP21)
summary(test.kp21)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2  78.25   39.12   2.197  0.142
## Residuals        17 302.75   17.81
samptemps.KP22 <- samptemps %>% filter(site == "Kebler Pass", year == "2022")
test.kp22 <- aov(avg.temp ~ as.factor(round), data = samptemps.KP22)
summary(test.kp22)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2   99.6   49.78   1.409  0.271
## Residuals        17  600.4   35.32
samptemps.OPH22 <- samptemps %>% filter(site == "Ohio Pass High", year == "2022")
test.oph22 <- aov(avg.temp ~ as.factor(round), data = samptemps.OPH22)
summary(test.oph22)
##                  Df Sum Sq Mean Sq F value Pr(>F)  
## as.factor(round)  2  66.87   33.44   4.459 0.0277 *
## Residuals        17 127.47    7.50                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
samptemps.OPL22 <- samptemps %>% filter(site == "Ohio Pass Low", year == "2022")
test.opl22 <- aov(avg.temp ~ as.factor(round), data = samptemps.OPL22)
summary(test.opl22)
##                  Df Sum Sq Mean Sq F value  Pr(>F)   
## as.factor(round)  2  162.2   81.11    6.06 0.00972 **
## Residuals        18  240.9   13.39                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
samptemps.RG21 <- samptemps %>% filter(site == "Rustlers Gulch", year == "2021")
test.rg21 <- aov(avg.temp ~ as.factor(round), data = samptemps.RG21)
summary(test.rg21)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2  23.24  11.621   1.563  0.237
## Residuals        18 133.82   7.434
samptemps.RG22 <- samptemps %>% filter(site == "Rustlers Gulch", year == "2022")
test.rg22 <- aov(avg.temp ~ as.factor(round), data = samptemps.RG22)
summary(test.rg22)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2  26.69  13.343   2.097  0.149
## Residuals        20 127.25   6.363
samptemps.SF21 <- samptemps %>% filter(site == "Stupid Falls", year == "2021")
test.sf21 <- aov(avg.temp ~ as.factor(round), data = samptemps.SF21)
summary(test.sf21)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2  14.04    7.02   0.435  0.654
## Residuals        17 274.08   16.12
samptemps.SF22 <- samptemps %>% filter(site == "Stupid Falls", year == "2022")
test.sf22 <- aov(avg.temp ~ as.factor(round), data = samptemps.SF22)
summary(test.sf22)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(round)  2  31.71   15.85   1.141  0.339
## Residuals        20 277.90   13.89

Round does significantly impact temperature, so it will make sense to include round as a main effect in our subsequent analyses.

4.2 Multinomial Choice

We use a mixed conditional logit model (mclogit) to determine whether pollinators’ foraging preferences are affected by temperature. In this model, we can imagine each interaction as a discrete choice made by a pollinator among the choice set of all flowers available at that site/date combination (defined by site.date). We include round as a random intercept and the identities of the insect and plant involved in the interaction as random slopes.

4.2.1 Functional Groups

First, we run this analysis with insects grouped by functional group. This comprises the largest amount of data.

mc.func <- mclogit(cbind(int, as.integer(date.site)) ~
                     poly(avg.temp, 2, raw = TRUE) +
                      as.factor(round),
                   random = list(~1 + avg.temp|gross_ID,
                                 ~1 + avg.temp|plant),
                   data = final.func)
## 
## Iteration 1 - deviance = 7609.395 - criterion = 0.9746249
## Iteration 2 - deviance = 7376.174 - criterion = 0.02211211
## Iteration 3 - deviance = 7343.175 - criterion = 0.00640909
## Iteration 4 - deviance = 7331.056 - criterion = 0.002133236
## Iteration 5 - deviance = 7327.29 - criterion = 0.0002740379
## Iteration 6 - deviance = 7326.449 - criterion = 6.167405e-06
## Iteration 7 - deviance = 7326.374 - criterion = 2.948911e-09
## converged
summary(mc.func)
## 
## Call:
## mclogit(formula = cbind(int, as.integer(date.site)) ~ poly(avg.temp, 
##     2, raw = TRUE) + as.factor(round), data = final.func, random = list(~1 + 
##     avg.temp | gross_ID, ~1 + avg.temp | plant))
## 
## Coefficents:
##                                 Estimate Std. Error z value Pr(>|z|)    
## poly(avg.temp, 2, raw = TRUE)1  0.224690   0.073511   3.057  0.00224 ** 
## poly(avg.temp, 2, raw = TRUE)2 -0.005419   0.001262  -4.294 1.76e-05 ***
## as.factor(round)2               0.106515   0.060862   1.750  0.08010 .  
## as.factor(round)3              -0.166317   0.070442  -2.361  0.01822 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Co-)Variances:
## Grouping level: gross_ID 
##          Estimate              Std.Err.           
## (Const.)  3.926835             2.935e+00          
## avg.temp -0.005289  0.004275   3.952e-03 5.306e-06
## 
## Grouping level: plant 
##          Estimate            Std.Err.       
## (Const.) 2.186126            1.03431        
## avg.temp 0.033240 0.004818   0.01576 0.00024
## 
## Approximate residual deviance: 7326 
## Number of Fisher scoring iterations:  7
## Number of observations
##   Groups by gross_ID: 15
##   Groups by plant: 64
##   Individual observations:  2585

4.2.2 Family

Next, we run this analysis with insects identified to the family level. This is a coarser-scale analysis than genus level, but there are a larger number of data points at the family level (because it includes insects for which we can identify to family but not to genus).

mc.family <- mclogit(cbind(int, as.integer(date.site)) ~
                     poly(avg.temp, 2, raw = TRUE) +
                       as.factor(round),
                   random = list(~1 + avg.temp|insect_family,
                                 ~1 + avg.temp|plant),
                   data = final.family)
## 
## Iteration 1 - deviance = 6469.196 - criterion = 0.9729606
## Iteration 2 - deviance = 6242.679 - criterion = 0.110585
## Iteration 3 - deviance = 6210.56 - criterion = 0.009696167
## Iteration 4 - deviance = 6198.641 - criterion = 0.002914789
## Iteration 5 - deviance = 6194.957 - criterion = 0.0003635539
## Iteration 6 - deviance = 6194.189 - criterion = 7.550973e-06
## Iteration 7 - deviance = 6194.122 - criterion = 2.866196e-09
## converged
summary(mc.family)
## 
## Call:
## mclogit(formula = cbind(int, as.integer(date.site)) ~ poly(avg.temp, 
##     2, raw = TRUE) + as.factor(round), data = final.family, random = list(~1 + 
##     avg.temp | insect_family, ~1 + avg.temp | plant))
## 
## Coefficents:
##                                 Estimate Std. Error z value Pr(>|z|)    
## poly(avg.temp, 2, raw = TRUE)1  0.170662   0.079078   2.158 0.030916 *  
## poly(avg.temp, 2, raw = TRUE)2 -0.004732   0.001390  -3.404 0.000664 ***
## as.factor(round)2               0.104609   0.069255   1.510 0.130917    
## as.factor(round)3              -0.106748   0.080308  -1.329 0.183770    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Co-)Variances:
## Grouping level: insect_family 
##          Estimate          Std.Err.           
## (Const.) 1.88751           1.6661393          
## avg.temp 0.03328 0.00509   0.0294493 0.0005205
## 
## Grouping level: plant 
##          Estimate            Std.Err.         
## (Const.) 1.971802            1.119256         
## avg.temp 0.038663 0.004551   0.021988 0.000432
## 
## Approximate residual deviance: 6194 
## Number of Fisher scoring iterations:  7
## Number of observations
##   Groups by insect_family: 37
##   Groups by plant: 64
##   Individual observations:  2124

We find that temperature has highly significant linear and quadratic effects.

4.2.3 Genus

Next, we run this analysis again but at the genus level. This is a finer-scale analysis that eliminates some uncertainty associated with the potential for variation among genera within the same family, but it loses some statistical power since we have roughly 500 fewer interactions that can be identified to the genus level.

mc.genus <- mclogit(cbind(int, as.integer(date.site)) ~
                     poly(avg.temp, 2, raw = TRUE) +
                      round,
                   random = list(~1 + avg.temp|insect_genus,
                                 ~1 + avg.temp|plant),
                   data = final.genus)
## 
## Iteration 1 - deviance = 5189.292 - criterion = 0.9178357
## Iteration 2 - deviance = 5086.866 - criterion = 0.02430653
## Iteration 3 - deviance = 5058.967 - criterion = 0.003680271
## Iteration 4 - deviance = 5029.556 - criterion = 0.001792912
## Iteration 5 - deviance = 5000.943 - criterion = 0.001973756
## Iteration 6 - deviance = 4974.1 - criterion = 0.001997291
## Iteration 7 - deviance = 4950.184 - criterion = 0.001860424
## Iteration 8 - deviance = 4929.955 - criterion = 0.001618803
## Iteration 9 - deviance = 4913.621 - criterion = 0.001343169
## Iteration 10 - deviance = 4900.911 - criterion = 0.001084788
## Iteration 11 - deviance = 4891.257 - criterion = 0.0008686005
## Iteration 12 - deviance = 4884.002 - criterion = 0.0007008919
## Iteration 13 - deviance = 4878.535 - criterion = 0.0005777871
## Iteration 14 - deviance = 4874.354 - criterion = 0.0004907416
## Iteration 15 - deviance = 4871.082 - criterion = 0.0004296845
## Iteration 16 - deviance = 4868.451 - criterion = 0.0003847052
## Iteration 17 - deviance = 4866.282 - criterion = 0.0003468096
## Iteration 18 - deviance = 4864.461 - criterion = 0.0003081801
## Iteration 19 - deviance = 4862.92 - criterion = 0.0002624472
## Iteration 20 - deviance = 4861.631 - criterion = 0.0002056703
## Iteration 21 - deviance = 4860.59 - criterion = 0.0001392979
## Iteration 22 - deviance = 4859.925 - criterion = 5.395226e-05
## Iteration 23 - deviance = 4859.545 - criterion = 2.051158e-05
## Iteration 24 - deviance = 4859.299 - criterion = 2.425215e-05
## Iteration 25 - deviance = 4859.164 - criterion = 3.435256e-05
summary(mc.genus)
## 
## Call:
## mclogit(formula = cbind(int, as.integer(date.site)) ~ poly(avg.temp, 
##     2, raw = TRUE) + round, data = final.genus, random = list(~1 + 
##     avg.temp | insect_genus, ~1 + avg.temp | plant))
## 
## Coefficents:
##                                 Estimate Std. Error z value Pr(>|z|)  
## poly(avg.temp, 2, raw = TRUE)1  0.162469   0.092899   1.749   0.0803 .
## poly(avg.temp, 2, raw = TRUE)2 -0.004019   0.001651  -2.435   0.0149 *
## round                          -0.029829   0.043986  -0.678   0.4977  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Co-)Variances:
## Grouping level: insect_genus 
##          Estimate                Std.Err.           
## (Const.)  3.026e-02              3.860e-06          
## avg.temp -1.238e-03  9.927e-05   1.582e-07 6.481e-09
## 
## Grouping level: plant 
##          Estimate            Std.Err.           
## (Const.)  4.41538            1.0734239          
## avg.temp -0.08288  0.00359   0.0201589 0.0003786
## 
## Approximate residual deviance: 4859 
## Number of Fisher scoring iterations:  25
## Number of observations
##   Groups by insect_genus: 52
##   Groups by plant: 64
##   Individual observations:  1627
## Note: Algorithm did not converge.

We still see a highly significant quadratic effect of temperature, but only marginal significance in the linear effect of temperature. This loss of significance seems likely to be related to the loss of statistical power in this genus-level analysis.

Since we found similar results at both the family and genus levels, we perform all subsequent analyses at the genus level to take advantage of the finer taxonomic resolution.

4.3 Relationship between Temperature and Pollinator Community

4.3.1 Abundance

We will check whether temperature affects the number of pollinators actively foraging during our sampling periods.

func.active <- final.func %>%
  mutate(round = as.factor(round)) %>%
  group_by(date, site, round, avg.temp, gross_ID) %>%
  summarise(num.active = sum(int))

total.func <- func.active %>%
  group_by(date, site, round) %>%
  summarise(total.active = sum(num.active))

func.subset <- func.active %>% 
  group_by(date, site, gross_ID) %>%
  tally() %>%
  filter(n == 3) %>%
  inner_join(func.active, by = c("date", "site", "gross_ID")) %>%
  dplyr::select(-n) %>%
  left_join(total.func, by = c("date", "site", "round"))

activity.func <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|gross_ID), 
                          data = func.subset, 
                         family = binomial)
summary(activity.func)
##  Family: binomial  ( logit )
## Formula:          
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,  
##     2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |  
##     gross_ID)
## Data: func.subset
## 
##      AIC      BIC   logLik deviance df.resid 
##   1888.6   1927.6   -934.3   1868.6      353 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev.  Corr  
##  date     (Intercept) 1.090e-01 3.302e-01       
##  site     (Intercept) 3.265e-09 5.714e-05       
##  gross_ID (Intercept) 4.904e+00 2.215e+00       
##           avg.temp    4.341e-03 6.589e-02 -0.95 
## Number of obs: 363, groups:  date, 46; site, 6; gross_ID, 7
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)                    -0.9604050  1.6569274  -0.580    0.562
## poly(avg.temp, 2, raw = TRUE)1 -0.0229519  0.1025394  -0.224    0.823
## poly(avg.temp, 2, raw = TRUE)2  0.0005544  0.0017317   0.320    0.749
## round2                         -0.0753297  0.0815248  -0.924    0.355
## round3                         -0.0773739  0.0887164  -0.872    0.383
actfunc.resid <- simulateResiduals(activity.func, n = 250)
plot(actfunc.resid)

We see no significance of any fixed effects. Diagnostic plots all look good.

family.active <- final.family %>%
  mutate(round = as.factor(round)) %>%
  group_by(date, site, round, avg.temp, insect_family) %>%
  summarise(num.active = sum(int))

total.family <- family.active %>%
  group_by(date, site, round) %>%
  summarise(total.active = sum(num.active))

family.subset <- family.active %>%
  group_by(date, site, insect_family) %>%
  tally() %>%
  filter(n == 3) %>%
  inner_join(family.active, by = c("date", "site", "insect_family")) %>%
  dplyr::select(-n) %>%
  left_join(total.family, by = c("date", "site", "round"))

activity.family <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|insect_family), 
                          data = family.subset, 
                         family = binomial)

summary(activity.family)
##  Family: binomial  ( logit )
## Formula:          
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,  
##     2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |  
##     insect_family)
## Data: family.subset
## 
##      AIC      BIC   logLik deviance df.resid 
##   1176.5   1211.6   -578.3   1156.5      236 
## 
## Random effects:
## 
## Conditional model:
##  Groups        Name        Variance  Std.Dev.  Corr  
##  date          (Intercept) 2.168e-01 4.657e-01       
##  site          (Intercept) 3.246e-09 5.698e-05       
##  insect_family (Intercept) 7.700e+00 2.775e+00       
##                avg.temp    1.149e-02 1.072e-01 -0.98 
## Number of obs: 246, groups:  date, 40; site, 6; insect_family, 14
## 
## Conditional model:
##                                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)                    -0.3495175  2.1520437  -0.162    0.871
## poly(avg.temp, 2, raw = TRUE)1 -0.0362095  0.1420572  -0.255    0.799
## poly(avg.temp, 2, raw = TRUE)2  0.0002365  0.0024072   0.098    0.922
## round2                         -0.1470876  0.1092397  -1.346    0.178
## round3                          0.0963853  0.1209730   0.797    0.426
actfam.resid <- simulateResiduals(activity.family, n = 250)
plot(actfam.resid)

No significance of any fixed effects at the family-level. KS test is showing significance consistent with overdispersion. We do know that our data are overdispersed by the model will not converge with family = betabinomial. Overdispersion should artificially lower statistical significance (Type I error), so since we don’t see any statistical significance, we think it’s fine.

genus.active <- final.genus %>%
  mutate(round = as.factor(round)) %>%
  group_by(date, site, round, avg.temp, insect_genus) %>%
  summarise(num.active = sum(int))

total.genus <- genus.active %>%
  group_by(date, site, round) %>%
  summarise(total.active = sum(num.active))

genus.subset <- genus.active %>%
  group_by(date, site, insect_genus) %>%
  tally() %>%
  filter(n == 3) %>%
  inner_join(genus.active, by = c("date", "site", "insect_genus")) %>%
  dplyr::select(-n) %>%
  left_join(total.genus, by = c("date", "site", "round"))

activity.genus <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|insect_genus), 
                          data = genus.subset, 
                         family = binomial)
summary(activity.genus)
##  Family: binomial  ( logit )
## Formula:          
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,  
##     2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |  
##     insect_genus)
## Data: genus.subset
## 
##      AIC      BIC   logLik deviance df.resid 
##    721.1    751.6   -350.5    701.1      146 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance Std.Dev. Corr  
##  date         (Intercept)  0.67368 0.8208         
##  site         (Intercept)  0.14004 0.3742         
##  insect_genus (Intercept) 19.01941 4.3611         
##               avg.temp     0.02255 0.1502   -0.97 
## Number of obs: 156, groups:  date, 33; site, 6; insect_genus, 12
## 
## Conditional model:
##                                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)                     3.229362   3.385593   0.954    0.340
## poly(avg.temp, 2, raw = TRUE)1 -0.304036   0.223456  -1.361    0.174
## poly(avg.temp, 2, raw = TRUE)2  0.005317   0.003841   1.384    0.166
## round2                         -0.049612   0.164785  -0.301    0.763
## round3                          0.230728   0.190392   1.212    0.226
actgen.resid <- simulateResiduals(activity.genus, n = 250)
plot(actgen.resid)

No significance of any fixed effects at the genus level. As with the family-level diagnostics, we see significance in the KS test consistent with overdispersion. The genus-level analysis will converge with a betabinomial distribution but lack of significance remains the same.

4.3.2 Composition

Next, we will use matrix regression to examine how dissimilarity in community composition relates to dissimilarity in temperature. For this analysis, we will include site/year and round as controlling factors that limit how permutations are conducted. We do not include round as a fixed effect because the MRM does not work with categorical variables.

plants.sum <-  plants %>%
  filter(plant != "") %>%
  mutate(site.date = paste(site, date, sep = ".")) %>%
  group_by(site.date, plant) %>%
  summarise(total_flowers = sum(total_flowers))

func.polmatrix <- final.func %>%
  group_by(concat, gross_ID) %>%
  summarise(num.active = sum(int)) %>%
  pivot_wider(names_from = gross_ID,
              values_from = num.active,
              values_fill = 0) %>%
  arrange(concat) %>%
  column_to_rownames(var = "concat")

func.polbeta <- as.matrix(vegdist(func.polmatrix, method = "horn"))

func.explanatory <- data.frame(concat = rownames(func.polmatrix)) %>%
  separate(concat, into = c("date", "site", "round"), sep = "[.]", remove = F) %>%
  mutate(date = as.Date(date),
         round = as.numeric(round),
         site.date = paste(site, date, sep = ".")) %>%
  left_join(samptemps, by = c("date", "site", "round")) %>%
  mutate(doy = strftime(date, format = "%j"),
         year = strftime(date, format = "%Y")) %>%
  arrange(concat)

func.tempdiff <- as.matrix(dist(func.explanatory$avg.temp))
rownames(func.tempdiff) <- rownames(func.polmatrix)
colnames(func.tempdiff) <- rownames(func.polmatrix)

func.doydiff <- as.matrix(dist(func.explanatory$doy))
rownames(func.doydiff) <- rownames(func.polmatrix)
colnames(func.doydiff) <- rownames(func.polmatrix)

func.rounddiff <- as.matrix(dist(func.explanatory$round))
rownames(func.rounddiff) <- rownames(func.polmatrix)
colnames(func.rounddiff) <- rownames(func.polmatrix)

func.plantmatrix <- func.explanatory %>%
  dplyr::select(site.date, concat) %>%
  left_join(plants.sum, by = "site.date") %>%
  dplyr::select(!site.date) %>%
  pivot_wider(names_from = plant,
              values_from = total_flowers,
              values_fill = 0) %>%
  arrange(concat) %>%
  column_to_rownames(var = "concat")

func.plantbeta <- as.matrix(vegdist(func.plantmatrix, method = "horn"))

## Permutation Controls
func.concatdf <- data.frame(row = 1:nrow(func.polmatrix),
                        concat = rownames(func.polmatrix)) %>%
  separate_wider_delim(cols = concat, names = c("date", "site", "round"), delim = ".") %>%
  mutate(year = strftime(date, format = "%Y"),
         site.year = paste(site, year, sep = "."))

func.control <- how(Within(type = "free"),
                     plots = Plots(as.factor(func.concatdf$site.year)),
                     blocks = as.factor(func.concatdf$round))

## MRM
set.seed(15723)
func.mrm <- suppressWarnings(MRM.base(as.dist(func.polbeta) ~ 
                  as.dist(func.tempdiff) + 
                  as.dist(func.plantbeta) + 
                  as.dist(func.doydiff),
                method = "logistic",
                nperm = 1000,
                permute.control = func.control))
func.mrm
## $coef
##                                    b  pval
## (Intercept)             -0.324750229 0.001
## as.dist.func.tempdiff.   0.017162610 0.096
## as.dist.func.plantbeta.  0.557472037 0.031
## as.dist.func.doydiff.    0.008174757 0.001
## 
## $dev
## resid.dev  resid.df  dev.pval 
##  7494.812 20096.000     1.000
family.polmatrix <- final.family %>%
  group_by(concat, insect_family) %>%
  summarise(num.active = sum(int)) %>%
  pivot_wider(names_from = insect_family,
              values_from = num.active,
              values_fill = 0) %>%
  arrange(concat) %>%
  column_to_rownames(var = "concat")

family.polbeta <- as.matrix(vegdist(family.polmatrix, method = "horn"))

family.explanatory <- data.frame(concat = rownames(family.polmatrix)) %>%
  separate(concat, into = c("date", "site", "round"), sep = "[.]", remove = F) %>%
  mutate(date = as.Date(date),
         round = as.numeric(round),
         site.date = paste(site, date, sep = ".")) %>%
  left_join(samptemps, by = c("date", "site", "round")) %>%
  mutate(doy = strftime(date, format = "%j"),
         year = strftime(date, format = "%Y")) %>%
  arrange(concat)

family.tempdiff <- as.matrix(dist(family.explanatory$avg.temp))
rownames(family.tempdiff) <- rownames(family.polmatrix)
colnames(family.tempdiff) <- rownames(family.polmatrix)

family.doydiff <- as.matrix(dist(family.explanatory$doy))
rownames(family.doydiff) <- rownames(family.polmatrix)
colnames(family.doydiff) <- rownames(family.polmatrix)

family.rounddiff <- as.matrix(dist(family.explanatory$round))
rownames(family.rounddiff) <- rownames(family.polmatrix)
colnames(family.rounddiff) <- rownames(family.polmatrix)

family.plantmatrix <- family.explanatory %>%
  dplyr::select(site.date, concat) %>%
  left_join(plants.sum, by = "site.date") %>%
  dplyr::select(!site.date) %>%
  pivot_wider(names_from = plant,
              values_from = total_flowers,
              values_fill = 0) %>%
  arrange(concat) %>%
  column_to_rownames(var = "concat")

family.plantbeta <- as.matrix(vegdist(family.plantmatrix, method = "horn"))

## Permutation Controls
family.concatdf <- data.frame(row = 1:nrow(family.polmatrix),
                        concat = rownames(family.polmatrix)) %>%
  separate_wider_delim(cols = concat, names = c("date", "site", "round"), delim = ".") %>%
  mutate(year = strftime(date, format = "%Y"),
         site.year = paste(site, year, sep = "."))

family.control <- how(Within(type = "free"),
                     plots = Plots(as.factor(family.concatdf$site.year)),
                     blocks = as.factor(family.concatdf$round))

## MRM
set.seed(71523)
family.mrm <- suppressWarnings(MRM.base(as.dist(family.polbeta) ~ 
                  as.dist(family.tempdiff) + 
                  as.dist(family.plantbeta) + 
                  as.dist(family.doydiff),
                method = "logistic",
                nperm = 1000,
                permute.control = family.control))
family.mrm
## $coef
##                                     b  pval
## (Intercept)               0.097492788 1.000
## as.dist.family.tempdiff.  0.013756718 0.275
## as.dist.family.plantbeta. 0.506663633 0.034
## as.dist.family.doydiff.   0.008344169 0.001
## 
## $dev
## resid.dev  resid.df  dev.pval 
##  9335.089 19896.000     1.000
genus.polmatrix <- final.genus %>%
  group_by(concat, insect_genus) %>%
  summarise(num.active = sum(int)) %>%
  pivot_wider(names_from = insect_genus,
              values_from = num.active,
              values_fill = 0) %>%
  arrange(concat) %>%
  column_to_rownames(var = "concat")

genus.polbeta <- as.matrix(vegdist(genus.polmatrix, method = "horn"))

genus.explanatory <- data.frame(concat = rownames(genus.polmatrix)) %>%
  separate(concat, into = c("date", "site", "round"), sep = "[.]", remove = F) %>%
  mutate(date = as.Date(date),
         round = as.numeric(round),
         site.date = paste(site, date, sep = ".")) %>%
  left_join(samptemps, by = c("date", "site", "round")) %>%
  mutate(doy = strftime(date, format = "%j"),
         year = strftime(date, format = "%Y")) %>%
  arrange(concat)

genus.tempdiff <- as.matrix(dist(genus.explanatory$avg.temp))
rownames(genus.tempdiff) <- rownames(genus.polmatrix)
colnames(genus.tempdiff) <- rownames(genus.polmatrix)

genus.doydiff <- as.matrix(dist(genus.explanatory$doy))
rownames(genus.doydiff) <- rownames(genus.polmatrix)
colnames(genus.doydiff) <- rownames(genus.polmatrix)

genus.rounddiff <- as.matrix(dist(genus.explanatory$round))
rownames(genus.rounddiff) <- rownames(genus.polmatrix)
colnames(genus.rounddiff) <- rownames(genus.polmatrix)

genus.plantmatrix <- genus.explanatory %>%
  dplyr::select(site.date, concat) %>%
  left_join(plants.sum, by = "site.date") %>%
  dplyr::select(!site.date) %>%
  pivot_wider(names_from = plant,
              values_from = total_flowers,
              values_fill = 0) %>%
  arrange(concat) %>%
  column_to_rownames(var = "concat")

genus.plantbeta <- as.matrix(vegdist(genus.plantmatrix, method = "horn"))

## Permutation Controls
genus.concatdf <- data.frame(row = 1:nrow(genus.polmatrix),
                        concat = rownames(genus.polmatrix)) %>%
  separate_wider_delim(cols = concat, names = c("date", "site", "round"), delim = ".") %>%
  mutate(year = strftime(date, format = "%Y"),
         site.year = paste(site, year, sep = "."))

genus.control <- how(Within(type = "free"),
                     plots = Plots(as.factor(genus.concatdf$site.year)),
                     blocks = as.factor(genus.concatdf$round))

## MRM
set.seed(62223)
genus.mrm <- suppressWarnings(MRM.base(as.dist(genus.polbeta) ~ 
                  as.dist(genus.tempdiff) + 
                  as.dist(genus.plantbeta) + 
                  as.dist(genus.doydiff),
                method = "logistic",
                nperm = 1000,
                permute.control = genus.control))
genus.mrm
## $coef
##                                     b  pval
## (Intercept)               0.203625449 1.000
## as.dist.genus.tempdiff.  -0.008774277 0.518
## as.dist.genus.plantbeta.  0.627047758 0.018
## as.dist.genus.doydiff.    0.015189146 0.001
## 
## $dev
## resid.dev  resid.df  dev.pval 
##  11623.12  18911.00      1.00

4.4 Relationship between Temperature and Plant Community

We will check whether temperature affects the attractiveness of plants during our sampling periods. We will run this at all three levels of pollinator taxonomy again, just in case there are any differences associated with the different subsets of data.

# Functional Group
func.plants <- final.func %>%
  mutate(round = as.factor(round)) %>%
  group_by(date, site, round, avg.temp, plant) %>%
  summarise(num.active = sum(int)) %>%
  filter(num.active > 0)

totalplants <- func.plants %>%
  group_by(date, site, round) %>%
  summarise(total.active = sum(num.active))

plants.subset <- func.plants %>% 
  group_by(date, site, plant) %>%
  tally() %>%
  filter(n == 3) %>%
  inner_join(func.plants, by = c("date", "site", "plant")) %>%
  dplyr::select(-n) %>%
  left_join(totalplants, by = c("date", "site", "round"))

activityplants <- glmmTMB(cbind(num.active, total.active-num.active) ~ poly(avg.temp, 2, raw = TRUE) + round + (1|date) + (1|site) + (1 + avg.temp|plant), 
                          data = plants.subset, 
                         family = binomial)
summary(activityplants)
##  Family: binomial  ( logit )
## Formula:          
## cbind(num.active, total.active - num.active) ~ poly(avg.temp,  
##     2, raw = TRUE) + round + (1 | date) + (1 | site) + (1 + avg.temp |  
##     plant)
## Data: plants.subset
## 
##      AIC      BIC   logLik deviance df.resid 
##   1544.2   1581.8   -762.1   1524.2      308 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev. Corr  
##  date   (Intercept) 1.022156 1.01102        
##  site   (Intercept) 0.267738 0.51743        
##  plant  (Intercept) 6.595667 2.56820        
##         avg.temp    0.005845 0.07645  -0.94 
## Number of obs: 318, groups:  date, 47; site, 6; plant, 30
## 
## Conditional model:
##                                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)                     1.779903   1.888670   0.942    0.346
## poly(avg.temp, 2, raw = TRUE)1 -0.176164   0.130978  -1.345    0.179
## poly(avg.temp, 2, raw = TRUE)2  0.003366   0.002356   1.429    0.153
## round2                         -0.004230   0.102757  -0.041    0.967
## round3                         -0.043848   0.116240  -0.377    0.706
plants.resid <- simulateResiduals(activityplants, n = 250)
plot(plants.resid)

5 Figures

5.1 Table 1

table1 <- final.genus %>%
  group_by(site, year, date, round) %>%
  tally() %>%
  group_by(site, year) %>%
  tally(name = "rounds.sampled") %>%
  right_join(final.genus, by = c("site", "year"), multiple = "all") %>%
  group_by(site, year, rounds.sampled) %>%
  summarise(days.sampled = length(unique(date)))
table1
## # A tibble: 10 × 4
## # Groups:   site, year [10]
##    site             year  rounds.sampled days.sampled
##    <chr>            <chr>          <int>        <int>
##  1 Brush Creek Hill 2021              16            6
##  2 Brush Creek Hill 2022              15            6
##  3 Kebler Pass      2021              20            7
##  4 Kebler Pass      2022              18            7
##  5 Ohio Pass High   2022              19            7
##  6 Ohio Pass Low    2022              20            8
##  7 Rustlers Gulch   2021              21            7
##  8 Rustlers Gulch   2022              23            8
##  9 Stupid Falls     2021              20            7
## 10 Stupid Falls     2022              23            8

5.2 Figure 1

Relationship between round and temperature.

fig1 <- ggplot(samptemps %>% mutate(site.year = paste(site, year, sep = ".")),
               aes(x = site.year, y = avg.temp, fill = as.factor(round))) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(color = "grey20",
              alpha = 0.25,
              show.legend = F,
             position = position_jitterdodge(dodge.width = 0.75,
                                             jitter.width = 0.2)) +
  scale_fill_manual(values = brewer.pal(n = 8, "Dark2")[c(3, 5, 6)]) +
  scale_color_manual(values = brewer.pal(n = 8, "Dark2")[c(3, 5, 6)]) +
  labs(y = "Average Temperature (C)",
       fill = "Round") +
  theme(axis.text.x=element_text(angle = 60, hjust = 1),
        text = element_text(size = 20),
        axis.title.x = element_blank())
fig1

# ggsave("round-temp.png", fig1, width = 10, height = 8, units = "in")

5.3 Figures 2 & 3

General relationship between temperature and interaction occurrence.

First let’s calculate some broad-scale trends in interaction occurrence. Specifically, we will find the top 10 most frequently observed pollinator genera and the top 10 most frequently visited floral partners.

top.pols <- final.genus %>%
  group_by(insect_genus) %>%
  summarise(total.int = sum(int)) %>%
  slice_max(order_by = total.int, n = 10)

top.plants <- final.genus %>%
  group_by(plant) %>%
  summarise(total.int = sum(int)) %>%
  slice_max(order_by = total.int, n = 10)

top.fams <- final.family %>%
  group_by(insect_family) %>%
  summarise(total.int = sum(int)) %>%
  slice_max(order_by = total.int, n = 10)

top.func <- final.func %>%
  group_by(gross_ID) %>%
  summarise(total.int = sum(int)) %>%
  slice_max(order_by = total.int, n = 6)

top.ints <- final.genus %>%
  filter(insect_genus %in% top.pols$insect_genus,
         plant %in% top.plants$plant) %>%
  group_by(insect_genus, plant) %>%
  summarise(total.int = sum(int)) %>%
  filter(total.int >= 5)

top.ints2 <- final.func %>%
  group_by(int.ID) %>%
  summarise(total.int = sum(int)) %>%
  slice_max(order_by = total.int, n = 50)

5.3.1 Functional Group Level

# Bombus
bombus.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Bombus", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
  filter(group == "2")

bombus.df[bombus.df$plant == "Potentilla pulcherrima",]$x[which.max(bombus.df[bombus.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 28.506
bombus.df[bombus.df$plant == "Heliomeris multiflora",]$x[which.max(bombus.df[bombus.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 35.918
bombus.df[bombus.df$plant == "Senecio crassulus",]$x[which.max(bombus.df[bombus.df$plant == "Senecio crassulus",]$predicted)]
## [1] 32.188
bombus.df[bombus.df$plant == "Erigeron speciosus",]$x[which.max(bombus.df[bombus.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 33.221
bombus.plot <- ggplot(bombus.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = brewer.pal(n = 5, "Dark2")[c(1, 2, 3, 5)]) +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  coord_cartesian(ylim = c(0, 0.004)) +
  theme(text = element_text(size = 20),
        legend.position = "none")
bombus.plot

ggsave("bombus-predcurves.png", bombus.plot, width = 6, height = 6, units = "in")

# Syrphidae
syrphid.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Syrphidae", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
  filter(group == "2")

syrphid.df[syrphid.df$plant == "Potentilla pulcherrima",]$x[which.max(syrphid.df[syrphid.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 26.543
syrphid.df[syrphid.df$plant == "Heliomeris multiflora",]$x[which.max(syrphid.df[syrphid.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 33.914
syrphid.df[syrphid.df$plant == "Senecio crassulus",]$x[which.max(syrphid.df[syrphid.df$plant == "Senecio crassulus",]$predicted)]
## [1] 30.014
syrphid.df[syrphid.df$plant == "Erigeron speciosus",]$x[which.max(syrphid.df[syrphid.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 30.965
syrphid.df[syrphid.df$plant == "Pseudocymopterus montanus",]$x[which.max(syrphid.df[syrphid.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 24.545
syrphid.plot <- ggplot(syrphid.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  coord_cartesian(ylim = c(0, 0.004)) +
  theme(text = element_text(size = 20),
        legend.position = "none")
syrphid.plot

ggsave("syrphid-predcurves.png", syrphid.plot, width = 6, height = 6, units = "in")

# SBSB
sbsb.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Small black solitary bee", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
  filter(group == "2")

sbsb.df[sbsb.df$plant == "Potentilla pulcherrima",]$x[which.max(sbsb.df[sbsb.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 18.806
sbsb.df[sbsb.df$plant == "Heliomeris multiflora",]$x[which.max(sbsb.df[sbsb.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 26
sbsb.df[sbsb.df$plant == "Senecio crassulus",]$x[which.max(sbsb.df[sbsb.df$plant == "Senecio crassulus",]$predicted)]
## [1] 22.238
sbsb.df[sbsb.df$plant == "Erigeron speciosus",]$x[which.max(sbsb.df[sbsb.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 23.198
sbsb.df[sbsb.df$plant == "Pseudocymopterus montanus",]$x[which.max(sbsb.df[sbsb.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 17.089
sbsb.plot <- ggplot(sbsb.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  coord_cartesian(ylim = c(0, 0.004)) +
  theme(text = element_text(size = 20),
        legend.position = "none")
sbsb.plot

ggsave("sbsb-predcurves.png", sbsb.plot, width = 6, height = 6, units = "in")

# Other Fly
fly.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Other fly", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
  filter(group == "2")

fly.df[fly.df$plant == "Potentilla pulcherrima",]$x[which.max(fly.df[fly.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 19.758
fly.df[fly.df$plant == "Heliomeris multiflora",]$x[which.max(fly.df[fly.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 27.084
fly.df[fly.df$plant == "Senecio crassulus",]$x[which.max(fly.df[fly.df$plant == "Senecio crassulus",]$predicted)]
## [1] 23.34
fly.df[fly.df$plant == "Erigeron speciosus",]$x[which.max(fly.df[fly.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 24.306
fly.df[fly.df$plant == "Pseudocymopterus montanus",]$x[which.max(fly.df[fly.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 17.855
fly.plot <- ggplot(fly.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  coord_cartesian(ylim = c(0, 0.004)) +
  theme(text = element_text(size = 20),
        legend.position = "none")
fly.plot

ggsave("fly-predcurves.png", fly.plot, width = 6, height = 6, units = "in")

# Beetles
beetles.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Beetles", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
  filter(group == "2")

beetles.df[beetles.df$plant == "Potentilla pulcherrima",]$x[which.max(beetles.df[beetles.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 22.525
beetles.df[beetles.df$plant == "Senecio crassulus",]$x[which.max(beetles.df[beetles.df$plant == "Senecio crassulus",]$predicted)]
## [1] 26
beetles.df[beetles.df$plant == "Erigeron speciosus",]$x[which.max(beetles.df[beetles.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 27.084
beetles.df[beetles.df$plant == "Pseudocymopterus montanus",]$x[which.max(beetles.df[beetles.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 20.424
beetles.plot <- ggplot(beetles.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = brewer.pal(n = 5, "Dark2")[c(1, 3, 4, 5)]) +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  coord_cartesian(ylim = c(0, 0.004)) +
  theme(text = element_text(size = 20),
        legend.position = "none")
beetles.plot

ggsave("beetles-predcurves.png", beetles.plot, width = 6, height = 6, units = "in")

# Hemipteran
hemip.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Heliomeris multiflora")), plant = "Heliomeris multiflora")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Senecio crassulus")), plant = "Senecio crassulus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Erigeron speciosus")), plant = "Erigeron speciosus")) %>%
  bind_rows(bind_cols(ggpredict(mc.func, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(gross_ID = "Hemipteran", plant = "Pseudocymopterus montanus")), plant = "Pseudocymopterus montanus")) %>%
  filter(group == "2")

hemip.df[hemip.df$plant == "Potentilla pulcherrima",]$x[which.max(hemip.df[hemip.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 25.806
hemip.df[hemip.df$plant == "Heliomeris multiflora",]$x[which.max(hemip.df[hemip.df$plant == "Heliomeris multiflora",]$predicted)]
## [1] 33.225
hemip.df[hemip.df$plant == "Senecio crassulus",]$x[which.max(hemip.df[hemip.df$plant == "Senecio crassulus",]$predicted)]
## [1] 29.356
hemip.df[hemip.df$plant == "Erigeron speciosus",]$x[which.max(hemip.df[hemip.df$plant == "Erigeron speciosus",]$predicted)]
## [1] 30.361
hemip.df[hemip.df$plant == "Pseudocymopterus montanus",]$x[which.max(hemip.df[hemip.df$plant == "Pseudocymopterus montanus",]$predicted)]
## [1] 23.919
hemip.plot <- ggplot(hemip.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  coord_cartesian(ylim = c(0, 0.004)) +
  theme(text = element_text(size = 20),
        legend.position = "none")
hemip.plot

ggsave("hemip-predcurves.png", hemip.plot, width = 6, height = 6, units = "in")

legend.plot <- ggplot(hemip.df, aes(x = x, y = predicted, color = plant, linetype = plant)) +
  geom_line(linewidth = 1) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability",
       color = "Floral Partner",
       linetype = "Floral Partner") +
  coord_cartesian(ylim = c(0, 0.011)) +
  theme(text = element_text(size = 18),
        legend.position = "bottom")  +
  guides(color = guide_legend(nrow = 2),
         linetype = guide_legend(nrow = 2))
ggsave("predcurves-legend.png", legend.plot, width = 10, height = 6, units = "in")

5.4 Figure 4

# SBSB
## Halictidae
halictidae.df <- data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.family, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(insect_family = "Halictidae", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  filter(group == "2")

fam.halic <- final.family %>%
  filter(insect_family == "Halictidae") %>%
  uncount(weights = int)

halictidae.df[halictidae.df$plant == "Potentilla pulcherrima",]$x[which.max(halictidae.df[halictidae.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 17.089
mcplot.halictidae <- ggplot(halictidae.df,
                            aes(x = x)) +
  geom_line(aes(y = predicted),
            color = "#7570B3",
            linewidth = 1.5) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "grey",
              color = "grey",
              alpha = 0.1,
              linetype = "dotted") +
  geom_vline(xintercept = 17.089,
             color = "black",
             linetype = "dashed") +
  geom_rug(data = fam.halic,
           aes(x = avg.temp),
           sides = "b") +
  coord_cartesian(xlim = c(15, 38),
                  ylim = c(0, 0.006)) +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability",
       color = "Plant",
       fill = "Plant") +
  theme(text = element_text(size = 20),
        legend.position = "bottom")
mcplot.halictidae

# ggsave("mcplot-halictidae.png", mcplot.halictidae, width = 6, height = 6, units = "in")

## Andrenidae
andrenidae.df <-   data.frame() %>%
  bind_rows(bind_cols(ggpredict(mc.family, terms = c("avg.temp [all]", "round [all]"), type = "fixed", condition = c(insect_family = "Andrenidae", plant = "Potentilla pulcherrima")), plant = "Potentilla pulcherrima")) %>%
  filter(group == "2")

fam.andren <- final.family %>%
  filter(insect_family == "Andrenidae") %>%
  uncount(weights = int)

andrenidae.df[andrenidae.df$plant == "Potentilla pulcherrima",]$x[which.max(andrenidae.df[andrenidae.df$plant == "Potentilla pulcherrima",]$predicted)]
## [1] 19.948
mcplot.andrenidae <- ggplot(andrenidae.df, aes(x = x, )) +
  geom_line(aes(y = predicted),
            color = "#7570B3",
            linewidth = 1.5) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "grey",
              color = "grey",
              alpha = 0.1,
              linetype = "dotted") +
  geom_vline(xintercept = 19.948,
             color = "black",
             linetype = "dashed") +
  geom_rug(data = fam.andren,
           aes(x = avg.temp),
           sides = "b") +
  coord_cartesian(xlim = c(15, 38),
                  ylim = c(0, 0.006)) +
  labs(x = "Temperature (C)",
       y = "Predicted Interaction Probability") +
  theme(text = element_text(size = 20))
mcplot.andrenidae

# ggsave("mcplot-andrenidae.png", mcplot.andrenidae, width = 6, height = 6, units = "in")

5.5 Figure 5

# Kebler Pass, 2021-07-21
kp.072121.bombus <- final.genus %>%
  filter(insect_genus == "Bombus",
         site == "Kebler Pass",
         date == "2021-07-21") %>%
  group_by(site, date, round, insect_genus) %>%
  summarise(genus.int = sum(int)) %>%
  left_join(final.genus, by = c("site", "date", "round", "insect_genus")) %>%
  group_by(round, insect_genus, plant, avg.temp, genus.int) %>%
  summarise(int = sum(int),
            rel.int = int/genus.int) %>%
  filter(rel.int > 0)

kp.072121.plant <- plants %>%
  group_by(date, site) %>%
  summarise(total_site = sum(total_flowers)) %>%
  left_join(plants, by = c("date", "site")) %>%
  filter(site == "Kebler Pass",
         date == "2021-07-21") %>%
  group_by(site, plant) %>%
  summarise(rel.abund = total_flowers/total_site) %>%
  mutate(plant = replace(plant, plant %in% c("Agoseris glauca", "Boechera stricta", "Collomia linearis", "Linum lewisii", "Helianthella quinquenervis", "Potentilla pulcherrima"), "0ther"))

kp.072121.plot <- ggplot(kp.072121.bombus,
                      aes(x = avg.temp, y = rel.int, fill = plant)) +
  geom_bar(position = "stack",
           stat = "identity",
           col = "grey50",
           width = 0.6) +
  scale_x_continuous(breaks=c(24.0, 30.7, 31.6)) +
  scale_fill_manual(values = brewer.pal(n = 9, "Blues")[c(3, 6, 8, 9)]) +
  labs(x = "Sampling Temperature (C)",
       y = "Interaction Occurrence",
       col = "Plant") +
  theme(text = element_text(size = 20),
        legend.position = "bottom",
        axis.title = element_blank(),
        legend.title = element_blank()) +
  coord_cartesian(xlim = c(22, 32))
kp.072121.plot

kp.072121.plantplot <- ggplot(kp.072121.plant,
                              aes(x = site, y = rel.abund, fill = plant)) +
  geom_bar(position = "stack",
           stat = "identity",
           width = 1.2) +
  scale_fill_manual(labels = c("Other", sort(unique(kp.072121.plant$plant))[2:5]),
                    values = c("lightgrey", "#C6DBEF", "#4292C6", "#08519C", "#08306B"))

# Rustlers Gulch, 2021-07-27
rg.072721.bombus <- final.genus %>%
  filter(insect_genus == "Bombus",
         site == "Rustlers Gulch",
         date == "2021-07-27") %>%
  group_by(site, date, round, insect_genus) %>%
  summarise(genus.int = sum(int)) %>%
  left_join(final.genus, by = c("site", "date", "round", "insect_genus")) %>%
  group_by(round, insect_genus, plant, avg.temp, genus.int) %>%
  summarise(int = sum(int),
            rel.int = int/genus.int) %>%
  filter(rel.int > 0)

rg.072721.plant <- plants %>%
  group_by(date, site) %>%
  summarise(total_site = sum(total_flowers)) %>%
  left_join(plants, by = c("date", "site")) %>%
  filter(site == "Rustlers Gulch",
         date == "2021-07-27") %>%
  group_by(site, plant) %>%
  summarise(rel.abund = total_flowers/total_site) %>%
  mutate(plant = replace(plant, plant %in% c("Artemisia dracunculus", "Erigeron elatior", "Potentilla fruticosa", "Pseudocymopterus montanus", "Senecio bigelovii", "Valeriana edulis"), "0ther"))

rg.072721.plot <- ggplot(rg.072721.bombus %>%
                          filter(insect_genus == "Bombus"),
                      aes(x = avg.temp, y = rel.int, fill = plant)) +
  geom_bar(position = "stack",
           stat = "identity",
           col = "grey50",
           width = 0.6) +
  scale_x_continuous(breaks=c(22.5, 25.1, 29.1)) +
  scale_fill_manual(values = brewer.pal(n = 9, "Blues")[c(1, 2, 5, 7)]) +
  labs(x = "Sampling Temperature (C)",
       y = "Interaction Occurrence",
       col = "Plant") +
  theme(text = element_text(size = 20),
        legend.position = "bottom",
        axis.title = element_blank(),
        legend.title = element_blank()) +
  coord_cartesian(xlim = c(22, 32))
rg.072721.plot

rg.072721.plantplot <- ggplot(rg.072721.plant,
                              aes(x = site, y = rel.abund, fill = plant)) +
  geom_bar(position = "stack",
           stat = "identity",
           width = 1.2) +
  scale_fill_manual(labels = c("Other", sort(unique(rg.072721.plant$plant))[2:6]),
                    values = c("lightgrey", "#F7FBFF", "#DEEBF7", "#C6DBEF", "#6BAED6", "#2171B5"))